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Carbon-based systems are prominent candidates for a solid-state spin-qubit due to weak spin-orbit 
and hyperfine interactions in combination with a low natural abundance of spin-carrying isotopes. 
We consider the effect of the hyperfine interaction on the coherence of an electron spin localized in 
a graphene quantum dot. It is known, that the hyperfine interaction in these systems is anisotropic 
promising interesting physics. We calculate the dynamics of an electron spin surrounded by a bath 
of nuclear spins in a non-Markovian approach using a generalized master equation. Considering 
a general form of the hyperfine interaction, we are able to extend the range of validity of our 
results to other systems beyond graphene. For large external magnetic fields, we find within Born 
approximation that the electron spin state is conserved up to small corrections, which oscillate with 
a frequency determined by the hyperfine interaction. The amplitude of these oscillations decays with 
a power law, where its initial value depends on the specific form of the anisotropy. Analyzing this 
in more detail, we identify two distinct classes of anisotropy, which can be both found in graphene 
depending on the orientation of the external magnetic field with respect to the carbon layer. 
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I. INTRODUCTION 

Quantum dots (QDs) in solid state nanostructures 
have attracted a lot of interest in recent years, partic- 
ularly, since localized spins are promising candidates'^ 
for qubits in spin-based quantum computers. Most of 
these devices are realized in semi-conducting compounds 
of III-V materials implying an environment of many nu- 
clear spins. In these structures, the hyperfine interaction 
(HI) between the nuclear spins and a central spin, for in- 
stance an electron or a hole spin, is a major source of both 
relaxation and decoherence of the central spinP^The dy- 
namics of these kind of spin systems has been extensively 
investigated both theoreticalljEH^I] anc j experimentally^ 
covering the control of the nuclear environmenpSlSU as 
well as the manipulation of the central spinP2H2U Since 
these investigations exhibit a decisive role of the HI, the 
use of materials with a low abundance of spin-carrying 
isotopes could offer a way to improve the properties of 
the qubits. 

Natural candidates for this are carbon and silicon 
based nanostructures, because the natural abundance of 
spin carrying 13 C and 29 Si is only about 1% and 5%, re- 
spectively. The experimental realization of carbon based 
QDs was achieved in recent years in graphen e " 51 ! as well 
as in carbonnanotubes^HSSl. i n silicon, the qubits can be 
fabricate d 56 * 57 ^ either using donor impurities or by confin- 
ing a single electron via electro-statical gates. However, 
a controlled localization of the donor impurities is still 
a challenging task and electro-statically confined Si QDs 
often involve nanostructures with other materials like Ge, 
which potentially introduce additional nuclear spins; see, 
for instance, Ref. [57] and references therein. 

Comparing the typical energy scale of the HI in 
different realizations of qubits reveals another advan- 
tage of carbon based systems. The hyperfine cou- 



pling constant! 58 ! 59 ! 61 ! Ai3c is significantly smaller than 
in GaAs, by about two orders of magnitude. Moreover, 
it is even less than the corresponding constant in Si based 
systems, both donor impurities^! and confined QDsJ^by 
approximately one order of magnitude. The s ame i s true 
for GaAs QDs using the spin of a heavy holeP^nSH 

In this article, we will investigate the properties of an 
electron localized in a QD fabricated from a graphene 
sheet, where the HI is anisotropic due to the p-type na- 
ture of the electrons in the conduction band. This type 
of interaction has not been studied extensively so far, be- 
cause the electronic properties of most qubits using the 
electron spin are primarily of the s-type and, hence, these 
systems are governed by an isotropic HI. Moreove r, t he 
few theoretical studies investigating anisotropic HP'S so 
far consider only a special case, where the longitudinal 
HI in the z direction deviates from the transverse in- 
teraction. Investigations of GaAs qubits using the spin 
of a heavy-hole analyzed a totally anisotropic HPS! or a 
HI, which is very weak in the transverse direction.^ Fi- 
nally, while most of the literature consi ders a n isotropic 
HI in Si systems, there are also studiegSSESJ taking the 
anisotropic HI into account. In contrast to these exam- 
ples, we will analyze a system with a HI, which is in all 
spatial directions on the same order of magnitude. Ex- 
perimentally, this physics is realized in a graphene QD, 
where the anisotropy can be varied by means of an ex- 
ternal magnetic field. 

This article is organized as follows: In Sec. 
most important results are summarized. In Sec. 
investigate the physical properties of our system and its 
Hamiltonian in detail. This part is followed by a brief 
recapitulation of the Nakajima-Zwanzig master equation 
in Sec. |IV[ which is used to investigate the dynamics 
of the electron spin in Sec. [V] Finally, we conclude in 
section Sec. IVH 
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II. KEY RESULTS 

We study the dynamics of an electron spin S in a 
graphene QD, where it is in contact with a bath of many 
nuclear spins located at sites via the HI as illus- 
trated in Fig. [l] Since there is no gap in the energy 
spectrum of graphene, the confinement of an electron is 
more challenging in contrast to electro-statically defined 
QDs in semi- conducting materials. One possibility to 
construct graphene nanostructures involves a chemical 
or mechanical treatment of graphene flakes, which has 
already been realized experiment all jBSHUl Furthermore, 
QDs can also be built by means of electrostatic potentials 
in presence of a finite energy gap, which can, for instance, 
be induced in single layer graphene by the substrate or in 
bilayer graphene by applying different potentials to the 
layers P2HLD 

We assume that the graphene flake is flat throughout 
the spatial extent of the QD, justifying our neglect of the 
influence of spin-orbit interaction on our problem. For 
simplicity, we consider further a rotational symmetric QD 
with the electron sitting in its ground state. Due to the 
confinement, the spatial probability distribution |</>(r)| 2 
of the electron is non-uniform and, thus, the HI becomes 
effectively site-dependent, since the major contribution 
to the coupling constant A arises from on-site terms Ak oc 
A- \4>(rk) | 2 - For clarity, we omit constant factors and refer 



to Eq. (151 for an exact definition of these couplings. 



The total action of the nuclear spins can be interpreted 
as a nuclear magnetic field 



J2 A 



(i) 



whose energy scale depends on the nuclear polarization 
p and the hyperfine coupling constant: (h z ) oc p ■ A. 
Additionally we allow for an external magnetic field B = 
(0, 0, B z ) giving rise to a Zeeman splitting b$ = frjsBz of 
the electronic energy levels, where 7s ~ 1.76 TO 11 T~ 1 s~ 1 
is the electron gyromagnetic ratio. This splitting is as- 
sumed to be much larger than the hyperfine energy: 
bg 3> A. In principle the magnetic field can be arbi- 
trary oriented with respect to the quantum dot plane. 
In this article, however, we restrict the orientation to be 
either perpendicular or parallel to this plane. Finally, 
we consider the system to be at zero temperature, which 
in particular implies that the thermal energy is much 
smaller than the spacing AE of the QD energy levels 
and the electron Zeeman energy: ksT <C bs, AE. 

This model is similar to previous studies of an elec- 
tron spin confined to a GaAs-QDp but there are, how- 
ever, graphene-specific characteristics, which lead to new 
physics. Since the natural abundance ri\ of spin-carrying 
isotopes is small for carbon, only N = nj ■ N tot of all 
atoms Ntot within the graphene dot carry spin, whereas 
all isotopes possess a spin in many semi-conducting ma- 
terials like GaAs, which is the most common material 
for building spin qubits. A comparison of the most im- 
portant properties of graphene and GaAs is given in Ta- 
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FIG. 1: (Color online) Schematic illustration of a two- 
dimensional quantum dot. An electron spin S is localized 
in a rotational symmetric QD, where it is in contact with 
bath of nuclear spins Due to the confinement, the spatial 
distribution of the electron in the ground state of the QD may 
be described by a Gaussian envelope function 0(r) given in 
Eq. (131, which in turn leads to a non-uniform HI between 
the electron spin and the nuclear spins. Finally, an external 
magnetic field is applied either perpendicular (B_l) or parallel 
(B||) with respect to the plane of the dot. 



ble[T] The HI coupling constant Ai3c estimate d 58 * 59 ! for 
ri/ = 1, is about two orders of magnitude smaller than 
the constant A(j a ^ s in GaAs lowering the nuclear mag- 
netic field by the same amount. For lower abundances 
the energy scale of this so-called Overhauser field is even 
further reduced by \h z ) oc AoaAs nj- Ai3 C . In the fol- 
lowing, we will explicitly distinguish between these two 
constants only if it is necessary to avoid confusion, while 
we use the general constant A apart from that. The com- 
parably small hyperfine energy in graphene extends the 
typical time scale thi oc h/A of this interaction signifi- 
cantly as compared to GaAs. Moreover, the use of much 
weaker external magnetic fields B becomes possible in a 
graphene QD, which can be quantified analyzing the ra- 
tio of the HI coupling constant and the electron Zeeman 
energy: 



A ~ niA I b s ■ 



(2) 



where we omitted factors of order O(l) for brevity and 
clarity. The exact definition is postponed to Eq. ( 32 ) in 
Sec. IV, where this ratio will serve as a small parame- 
ter justifying our perturbative treatment of the HI. For 
a natural 13 C-concentration of nj = 0.01, the Zeeman 
energy due to a comparably weak external magnetic field 
|B| = 2.6 mT still exceeds the hyperfine energy by more 
than two orders of magnitude: A « 10 -2 . 



Besides this experimentally motivated advantages, the 
study of a graphene QD is of general interest, since it 
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■ 2.6 • 10" 


thi oc h/A 




1 




100 



type of HI 

(X X , Xy, X Z ) 



(1,1,1) ±: (-1/2,-1/2,1) 
||: (1,-1/2,-1/2) 



TABLE I: Comparison of the most important parameters of 
GaAs and graphene. The total number of nuclei N to t is esti- 
mated for a QD of typical size R = 50 nm. While all nuclei 
in GaAs carry spin, the abundance ni of 13 C can in principle 
be modified, where the natural abundance is only nj = 0.01. 
The HI constant A in GaAs is about two orders of magnitude 
higher than in graphene (with ni — 1) demanding lower ex- 
ternal magnetic fields B z and leading to a prolonged typical 
hyperfine time scale thi- The specific values of the anisotropy 
constants Xi in graphene depends on the orientation of the 
external magnetic field, while in GaAs the HI is isotropic for 
electrons. 



demands the investigation of an anisotropic HI, which 
arises from the p-type nature of the electrons: 



Hhi = X z h z S z + X x h x S x 



XyhySy 



(3) 



This anisotropy is expressed in terms of three coupling 
constants X x ,X y , and \ z , where the isotropic limit is 
given by X x — X y — X z as it is, for instance, found 
for an electron in GaAs. In graphene, this anisotropy 
is related to the geometry of the carbon plane and, 
hence, arises between the in-plane (x, y) and out-of- 
plane (z) components^. Defining the quantization axis 
z by the external magnetic field B, the specific form 
of the anisotropy in graphene for different orientations 
is obtained by a geometric projection: {x, y, z} — ¥ 
{x,y, z}. This generally allows for a wide range of pos- 
sible anisotropies. In this article, we focus on two spe- 
cial cases, where a magnetic field perpendicular to the 
graphene plane gives rise to — 2X X = —2X y = X z = 1, 
while an in plane magnetic field generates another form 
of anisotropy: X x — —2X y = —2X Z = 1. 

Due to the large Zeeman splitting, real spin flips be- 
tween the electron and an arbitrary nucleus are forbidden 
by energy conservation, since the resulting energy dif- 
ference cannot be compensated by the hyperfine energy. 
Thus, only virtual spin flip-flops, where the electron spin 
is flipped back and forth, are allowed. If the characteris- 
tic time scale Tdd for dipolar spin interactions within the 
nuclear bath is much longer than the typical time scale 
thi of the HI, these virtual processes are the only source 
of change in the configuration of the nuclear spins. As a 



consequence of this, there is no randomizing effect within 
the nuclear bath itself and, hence, the system is in a non- 
Markovian regime rather than in a Markovian regime, 
where a stochastic nuclear magnetic field h would lead 
to a fast decay of the electron spin. However, even in a 
non-Markovian regime arbitrary states of the nuclear spin 
bath will in general lead to a fast decrease of the electron 
spin amplitude^!. Thus, the nuclear spins have to be pre- 
pared in a so-called narrowed state, which is an eigenstate 
of the nuclear magnetic field: h z \n) — \h z / \ n )- Such 
a narrowing can in principle be obta ined by a measure- 
ment of the nuclear magnetic fieldi2ttl2EU or by polariz- 
ing the bath to a high degree p ~ 1, which is, however, 
hardly feasible in experiments. Hence, the experimen- 
tally most promisi ng approaches are realized by pump- 
ing scheme j. 24 l 25 l 25 l t33l ; where the HI itself is employed to 
achieve such a state-narrowing. 

Assuming a narrowed nuclear spin state, the non- 
Markovian dynamics of the electron spin can be ana- 
lyzed for arbitrary anisotropy as defined in Eq. ([3|. The 
calculations are carried ou t by the use of the so-called 
Nakajima-Zwanzig equatio n 1 9 * 22 * 72 * 73 ^, which provides an 
integro-differential equation for the electron subsystem 
only. This equation, however, can not be solved ex- 
actly, but demands for further simplifications, which are 
obtained by a second order expansion in the HI. This 
truncation of higher order terms, however, sets an upper 
bound to the time regime t -C A -1 t#/, in which our 
predictions are valid. For a detailed description of this 
expansion, we refer to Section [TV| and turn to the results 
of these calculations. 

The spin dynamics will be described by expectation 
values of the transverse and longitudinal components, 
(S + )(t) and (S z )(t), respectively, which allow the inter- 
pretation of our results in terms of Bloch-like equations 
of motion. 

The expectation value of the transverse electron spin 
component exhibits two oscillating contributions, which 
differ in both their amplitudes and frequencies: 



(S+)(t) = a osc {t) 
The first term is given by 

( S +)o 



pow 
dec 



(*)• 



-M 



1 + (A 2 + A 2 )/4 • 6 



exp 



1 T* 



(4) 



(5) 



which describes a simple precession with a frequency de- 
termined by the Zeeman energy to n of the effective mag- 
netic field. This field consists of the external magnetic 
field B z and the nuclear magnetic field h z pointing in the 
same direction. The second contribution oscillates with 
a much smaller frequency t h ^, where the HI timescale 
thi = 2NH/\X z \n I A is depending on both the HI energy 
scale A and the absolute value of the longitudinal cou- 
pling constant |A Z |. Moreover, this second term is much 
smaller of order 5 = A 2 /N -c A and exhibits a power- 
law decay of its amplitude. The reason for the smallness 
of this contribution is the large Zeeman-splitting, which 
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makes the virtual flip-flop processes very unlikely due to 

the enormous energy cost of the electron spin flipping. 

For times t < t^i only a numerical solution for <r p ^(t) 

1 



was within reach, whereas for times thi t A 
this function is asymptotically described by: 



thi 



0(5), given by the initial value (S z ) ~ \Sz) - The sec- 
ond term behaves similar to the decaying part of the 
transverse spin component o^°™(£) except for an addi- 
tional overlying oscillation with frequency ui n . More im- 
portantly the constants 



.THI I T j ■ , t 

o t/sm 

t 



thi 



iT^cos(- 



t 

thi 



(6) 



The constants U and V are determined by the initial 
value (S+) Q = (S x ) + i(Sy) i t ne polarization p of the 
nuclear bath and the couplings in the transverse plane, 
X x and A y , respectively: 



and 



U = \l(S X ) Q +i\ x (S y ) 



v = p-(xl(s x ) Q + i\l(s y ) ) 



(7) 



(8) 



These results are true for A z > 0, from which one obtains 
the time evolution of the electron spin for A z < by ap- 
plying the following relation: (— |A z |,o;) O (|A 2 |,— ui). 
Thus a negative coupling in longitudinal direction can 
be treated as an inversion of the effective nuclear mag- 
netic held. From Eqs. (J7]) and ([8]), we see, that there 
are two distinct classes of anisotropy, which differ in the 
amplitude of the decaying contribution er^f (i) . For the 
first class, defined by X x — X y , the constants U and V 
are both proportional to the initial value of the trans- 
verse spin \S+) = \S X ) + i(<Sj/) - The reason for this 
is the preserved rotational symmetry in the transverse 
plane, which is for instance realized in a GaAs QD with 
X x = X y = X z = 1 and in a graphene QD subjected to a 
perpendicular magnetic field with — 2X X = —2X y = X z — 
1. Due to this symmetry, there is no precisely defined 
direction within this plane distinguishing between (S x ) 
and (Sy} . As a consequence of this, only their super- 
position (S+) = (Sx) + i(^y)o ma tt ers f° r the initial 
state preparation. 

This is, however, not the case for the second class of 
anisotropy defined by X x =/= X y , which is for example 
found in a graphene QD with in-plane magnetic field, 
where the couplings are given by X x = —2X y = —2X Z = 1. 
Now, the constants U and V describe a weighted mixture 
of (S X ) Q and (S y ) Q , which is caused by the broken rota- 
tional symmetry. This fact, in turn, allows to precisely 
define the x- and y-direction, with respect to which the 
transverse electron spin component can be prepared. A 
more detailed discussion of this is given at the end of 
section V.B. 

In contrast to this, the longitudinal spin component 
does not show qualitative differences between the two 
classes. Its time evolution consists of two parts: 



(S z )(t) = (S z )+*? e :(t).exp 



(9) 



This expectation value is dominated by a constant con- 
tribution \Sz), which is, up to small correction of order 



and 



U 



V 



X X Xy 



X 2 ^. + X y 



X X X y 



XI 



xi 



A~ + A~ 



A 2 . + A2 



(10) 



(11) 



show only a quantitative dependence on the anisotropy, 
while they are qualitatively unchanged for arbitrary 
anisotropy. Thus, the longitudinal spin component shows 
the same qualitative behavior for any reasonable values 
of X x and X y . 

Finally, we focus on the implications of our general 
results on the two graphene models, where the external 
magnetic field is either perpendicular (_L) or parallel (||) 
with respect to the carbon layer. While the former case 
shows the same qualitative behavior as an isotropic sys- 
tem, one faces a mixing of the initial transverse ampli- 
tudes due to the broken rotational symmetry for the later 
case. Moreover, as we have shown above, the time regime 
for which we find only a small partial decay is determined 
by the time scale of the HI. To be more precise, most of 
the electron spin is preserved for times t <C A _1 r^f/. 
Since this timescale depends on the absolute value of the 
coupling constant |A Z | in the longitudinal direction, one 
observes a different behavior for perpendicular and par- 
allel magnetic field, which is quantitatively described by 
the following ratio: 



'HI 



•hi 



|4| 



(12) 



Thus, the time regime, on which only a small partial 
power law decay of the spin amplitude and, vice versa, a 
mostly preserved electron spin is to be expected, is pro- 
longed by a factor of 2 for a parallel magnetic compared 
to a perpendicular oriented field. Therefore, we conclude, 
that a parallel orientation is preferable, although one has 
to prepare the transverse electron spin more carefully at 
the beginning of an experiment. 



III. MODEL AND HAMILTONIAN 

A. Physical properties of the model 

In this section, we discuss some of the physical prop- 
erties of the model introduced at the beginning of Sec. 
[TT| in more detail. For brevity and clarity, we measure all 
spins in units of H and adjust the coupling constants to 
give the right dimensions. 
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The quantum dot hosting the electron is defined by 
electro-statical gates on top of a graphene sheet, where 
we assume a rotational symmetric QD for simplicity. Due 
to the confinement, the QD has a discrete spectrum of 
bound stat es, with an energy splitting between different 
states^EiHlll jf t ne temperature is small compared to 
the level-spacing AE of these bound state energies, the 
electron will occupy the ground state, which we describe 
by an envelope-function of the form: 



4>{v) = 



1 



/ttR 



exp 



1 / r \ 2 

2 \RJ 



(13) 



where r = |r| is the absolute value of the electron po- 
sition. This wave-function is a Gaussian with a Bohr 
radius R, which we will regard as the radius of our QD. 
Note that the envelope function is not the exact electron 
wave-function, but should give a good approximation to 
the precise solution. This can be seen, for instance, in 
graphene QDs based on semi-conducting armchair nano- 
ribbons-^1. The most important aspects, which are cap- 
tured by this specific choice are the absence of nodes in 
the ground-state, a peak of the wave-function in the cen- 
ter as well as a strong decay inside the barriers. A QD 
with radius R includes 



(14) 



nuclei in total, where Vq = \/3a 2 /A is the (two- 
dimensional) unit cell volume containing one nucleus. 
For a typical graphene QD of size R — 50 nm and the 
graphene lattice constant given by a — 2A6A, there are 
Ntot ~ 10 5 nuclei within the dot. 

Naturally only a fraction nj ps 0.01 of these nuclei 
is of the spin-carrying species 13 C, while the remainder 
consists of spinless 12 C. However, it should in princi- 
ple be possible to modify this natural abundance by ei- 
ther isotopic purification or 13 C enr ichment , as it was 
already done for other carbon systems ^ 2 * 79 * 80 !. This leads 
to a fraction nj = N/N to t of 13 C-atoms carrying spin 
/ = 1/2, where < nj < 1. Especially the lower bound, 
where the number of nuclei becomes too small to treat 
them as a bath, is an interesting limit, which is, however, 
beyond the scope of this article. Besides the 13 C spins, we 
neglect all other possible sources of nuclear spins stem- 
ming from the gate- and substrate-materials, which is 
justified by the strong confinement of the electrons into 
the graphene plane. Altogether, we consider N = niN tot 
nuclear spins 1=1/2 within the QD radius R forming a 
spin-bath with an arbitrary polarization p, where p = 
for an unpolarized bath and p = 1 for full polarization. 

We first turn to the interaction between a single nu- 
clear spin 1^ = (Ik.x, Ik.y, Ik.z), located at an arbitrary 
position rfc within the QD, and the electron spin S at 
position r. In gene ral, the HI consists of three ma- 
jor contributions^^ the Fermi contact interaction, the 
anisotropic HI, and the coupling of the electron orbital 
angular momentum to the nuclear spins. For a carbon 
nanotube one would have to take into account all three 



terms due to the curvature dependent mixture of s- and 
p-type wave-functions, while for a perfectly flat sheet of 
graphene only the anisotropic HI is relevaniP^, where the 
anisotropy can be formulated in terms of coupling con- 
stants Xi and the strength of the HI is given by another 
constant A. This anisotropy arises from the p-type na- 
ture of the electron wave function and is, thus, related to 
the geometry of the graphene plane, where the in-plane 
and out of plane components deviate. Thus, the defi- 
nition of the quantization axis via an external magnetic 
field allows us to create different forms of anisotropy as 
stated in Sec. |TTJ 

Furthermore, the HI also depends on the position of 
the nucleus due to the non-uniform spatial distribution of 
the electron within the QD, since the HI is dominated by 
on-site contributions, where the electron and the nucleus 
are located at the same site. Thus, the interaction Ak at 
site k is determined by the hyperfine coupling constant 
A and the electron envelope function at r^: 



A k =AVo\<f>(T k )\' 



(15) 



where Vq is the volume of a primitive unit cell. For a fully 
13 C enriched system nj = 1 the hyperfine constant ia^S 
A = 0.6 neV implying that a single scattering process 
is on the order of A/N tot = njA/N PI If we relate the 
radial position < R of an arbitrary nuclear spin to the 
number of nuclear spins k within this radius by 



k 

N' 



(16) 



and use Eq. (13 1, we can rewrite the expression of the 
HI: 



njA 



exp 



k 
"N 



(17) 



for a system of N = njN to t nuclear spins. 

One important degree of freedom of an electron in 
graphene has so far not been taken into account. The 
model used here, does not include the valley degree of 
freedom. The valleys correspond to the two inequivalent 
Dirac cones of the hexagonal Brillouin zone and must be 
considered in scattering processes involving a short-range 
potential. As a consequence of this, one has to check, if 



the 



dependency of the anisotropic HI is sufficient 



to generate valley scattering. Therefore, we additionally 
investigate the momentum dependence of the coupling 
constants A = A (Q), where Q is the difference between 
incoming and outgoing momentum. Assuming, that the 
most important contributions arise from on-site terms^, 
we estimate, that the valley mixing Q = AK = K K' 
and valley conserving Q = processes are of comparable 
strength. 

This result reflects the short-range nature of the 
anisotropic HI leading to a broad distribution in mo- 
mentum space. Hence, the valley scattering is relevant, 
unless it is forbidden due to energy conservation, which 
demands that the HI compensates the energy splitting 
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of the two valleyfP. From Eq. (HtI, we know, that 



the energy scale of scattering processes between the elec- 
tron and a single nuclear spin is on average of order 
€hi <x~ njA/N 10~ 12 eV . In typical graphene QDs, 
we expect a non-zero valley splitting Ak,k' because of 
several mechanisms. First, every experimental setup will 
have deviations from an "ideal" setup, e.g. roughness 
of the boundaries or adatoms, which couple the valleys. 
But even in an ideal experiment, one can expect a valley- 
splitting due to the presence of an external magnetic- 
field in combination with a finite mass term induced 
by the substrate^!. Furthermore, if we assume a QD 
made of a semi-conducting graphene ribbon with arm- 
chair boundaries, the valley degeneracy is also lifted by 
h ^ F ~ lOmeV, where w = 50 nm is a typical width 
of a ribbon and vf — 10 6 to/s is the Fermi- velocity in 
graphene. 

Finally, let us briefly comment on the dynamics within 
the nuclear spin bath. This internal dynamics would be 
induced by nuclear dipolar interaction, which is, how ever, 
suppressed in QDs due to the "frozen core " effecPJSSl 
and, hence, negligible. As a consequence of this, we as- 
sume that the dynamics of the bath can be neglected on 
time scales t ~ tjji relevant for the electron spin evolu- 
tion. 



B. Hamiltonian 

We start this paragraph by formalizing the most im- 
portant facts of the two previous sections, where we have 
discussed the physical properties of our model in detail. 
Note that we replace classical variables S, I, h by the re- 
spective operators S,I, h henceforth. In the following 
section, we will consider an arbitrary anisotropy, which 
allows us to draw conclusions for a graphene QD sub- 
jected to either a perpendicular (_L) or a parallel (||) ori- 
ented external magnetic field. 

In any case, the respective field gives rise to a large 
Zeeman-splitting bs — fojsB z of the electron states, 
which in turn forbids real spin flip processes due to energy 
conservation. As a consequence of this, the Hamiltonian 
can be split in an unperturbed part H consisting of all 
Zeeman terms and a perturbative part Hy containing the 
virtual flip-flop processes of the HI. 

If we choose the quantization axis along the magnetic 
field, the unperturbed Hamiltonian reads: 



H = bil z + b s S z + X z h z S z 



(18) 



where the total z-component of the nuclear spins is given 
by I z = ^2 k Ik, z - Transforming to a rotating reference 
frame^, one can eliminate the nuclear Zeeman-term bji z 
in Eq. (18), which allows to write the unperturbed 



Hamiltonian Hq in compact form: 
H = QjS z , 



where the Zeeman-energy operator is given by 

u = b s - h + X z h z . (20) 

The perturbative part Hy of the Hamiltonian can be 
written as 



Hy — X X h X S X + XyhySy 

' g+S- + g-S. 



+ > 



(21) 



where we introduced raising and lowering operators S± — 
S x ± iS y and generalized nuclear operators g± , which are 
given by 



9± = 



(X x ± X y ) h + + (X x =F A y ) h- 
A±/i+ + Aip/i_ 



(22) 



Note, that the bare nuclear magnetic field operators h± — 
h x ± ih y can also be considered in terms of single nuclear 
spin raising and lowering operators located at sites k: 



From Eq. (22), it becomes obvious, that there are two 
classes of anisotropy with = A y and X x =^ X y , re- 
spectively. In the former case, the generalized operators 
g± oc h± are just multiples of the bare nuclear magnetic 
field operators, whereas in the later case these operators 
become linear combinations of both raising and lowering 
operators h±, giving rise to more flip-flop processes. 



C. Initial condition 

The time evolution of the combined system, consisting 
of the electron spin and N nuclear spins, is given by the 
action of the total Hamiltonian H = Hq + Hy on the 
initial state of the system at time t = 0, which is assumed 
to be a product state of the form: 



|^(0)) = |n)®|x) 



(23) 



The z-projection of the electron spin is either parallel 
(x = t) or anti-parallel (J,) to the magnetic field, while 
the nuclear spin bath is prepared in a so-called narrowed 
state at the beginning of an experi ment in order to pre- 
vent a fast decay of the electron spirP^. As discussed in 
detail in Sec. |Hj this narrowing can be created by mea- 
suring the nuclear spin system into an eigenstate of the 
z-component of the nuclear magnetic field h z : 



h z \n) = (h z ) |n) 



(24) 



(19) 



In general, this narrowed state is a superposition of many 
degenerate eigenstates \rij) having all the same polariza- 
tion p: 

I") = a i ' ^ = (^ z )n ' ( 25 ) 

3=1 
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Symbol 


Definition 


Symbol 


Definition 


SHI 


\\ z \mA/2N 


A 


niA 1 2oj n 


TH1 


2Nh / \\g\niA 


A 


A • (Xl + X 2 y ) / 2\X Z \ 


c± 


[lTP]/2 


<5 


A 2 / TV 


N 


N(\l + A 2 )/2A 2 


8 


S-(Xl + X 2 y )/2 



TABLE II: Important symbols used in the text. In the main 
part of this article, all energies (times) are measured in units 
of en i ijHi), which is the typical energy (time) scale of a sin- 
gle HI process. The parameters c± arise from expectation val- 
ues of the nuclear bath operators h± for spin 1/2 nuclei, which 
are characterized by their polarization p. The effective num- 
ber N is introduced to rewrite the denominator given in Eq. 
(61 1 in a generalized form for arbitrary anisotropy (X x , X y ,X z ). 
The total Zeeman energy is given by to n = bs — bj + X z (h z ) , 
where 6s (6/) is the electron (nuclear) Zeeman energy and 
( *) = P I n i A. For large external magnetic fields, the small 
parameter A 1 determines the perturbative regime, while 
the second small parameter S quantifies the non-Markovian 
corrections. The shorthands A and 5 allow a compact nota- 
tion of our results. 



The corresponding eigenvalue of the nuclear field opera- 
tor h z is given by^l 



(hz) n =pIniA. 



(26) 



This result is obtained by converting sums into integrals 
using a continuum limit, which is valid for times t <C 
y/N/2 thi, where the corrections to the integrals are 
still small. 

In the following, we introduce dimensionless quantities 
by measuring energies in units of ejjj = \X z \njA/2N, 
where N = njN to t is the actual number of nuclear spins 
within the dot. This typical energy scale for a single HI 
process also defines the corresponding time scale thi = 



£Hi/h = 2NH/\X z \niA, which serves as a measure for 
times henceforth. For clarity, we give a summarizing list 
of the most important symbols in Tab. [TXJ. 



IV. METHOD 

A. Nakajima-Zwanzig-equation 

A common tool to investigate systems like the one con- 
sidered in this article is the density matrix formalism. 
The time evolution of the total system is described by the 
von-Neumann equation p(t) = —\[H,p(t)] = —iCp{t), 
where the Liouville-(super-)operator C = Co + Cy also 
splits in a perturbed and an unperturbed part. The ini- 
tial density matrix is determined by the initial state of 
our system given in Eq. Q: p(0) = \i/)(Q)) {ip{0)\ = 
pi <8> ps (0). In order to find the dynamics of the electron 
subsystem only, one can rewrite this differential eq uation 
in terms of the Nakajima-Zwanzig master equatio n 1 72 ! 73 !, 
where a projection super-operator VO = pi-TiiO is used 
to perform a partial trace over the nuclear subsystem. 
This choice of the projector preserves both the initial 
condition 7^/5 (0) = p(Q) and the electron spin expec- 
tation valuesP^I (5/3) (t) = Tr[SpVp(t)\ = Tr[Spp(t)}, 
(3 — z,±. For later convenience, we also define the com- 
plement projector Q by V + Q = 1. Finally, the equation 
of motion can be written as 



i>s (t) = -iQp s (t)-i dt' S s (* - Ps (0 , (27) 



where the new Liouville operator CqOs = [oJ n ■ S z , 0$]- 
originates from the action of the Zeeman energy opera- 
tor uj on the initial nuclear state pr. ui n = T¥j (Cj ■ pi). 
While this first part describes the precession of the elec- 
tron spin in the effective magnetic field, the second term 
describes the flip flop processes arising from the HI. In 
order to dissolve the convolution between the self-energy 
Eg (t) and the spin density matrix, we perform a Laplace 
transformation to obtain a much simpler algebraic form: 



sps (s) - ps (0) = -i Co ps (s) - i (s) Ps (s) , 



(28) 



Ss(s) = -iTrj 



oo 



. 1=0 



1 

s + iQCo 



-ifiA 



1 

a + iQCo 



CvPi 



(29) 



3=1 



The series for the self-energy Eg (s) in powers of the 
interaction Liouvillian Cy is obtained by inserting the 
Laplace-transformed Schwinger-Dyson identity. For the 



explicit calculation of these self-energy terms, it is con- 
venient to choose a set of 4 basis vectors given by u^j. = 
~ (<T ± &g) and u ± = \ (& x ± & y ), where {oi} l=x y z are 



the 2 x 2-Pauli-matrices and <jq = 1 is the identity ma- 
trix. In this basis, the electron spin density matrix forms 
a 4-component vector 



second order, the self-energy takes the specific form: 



Ps = PT""T + PlH + (S+)u- + (S-)u+ 
= (P1-.PI>( S +)>( S -)) T ■ 



(30) 



The unperturbed Liouvillian is given by a diagonal ma- 
trix £o = 5 • diag (£_,—£_, — £+, £+), where the new 
operators 



c±6i = [w,6/] : 



(31) 



are defined by their action on an arbitrary nuclear spin 
operator O/. The perturbative Liouvillian Cy has a 4 x 4 
off-diagonal form containing the generalized nuclear mag- 
netic field operators g± . For further details of this calcu- 
lation, we refer to Appendix A of Ref. [9] and note that 
the bare nuclear magnetic field operators h± there cor- 
respond to the generalized operators g± of our model. 
In general, all contributions £( 2 -? +1 ) (s) to the self-energy 
containing odd powers of Cy vanish , since only virtual 
flip-flop processes are allowed by energy-conservation. 
Each even summand E^ 2 -? -1-2 ) (s) is associated with a small 
paramete^ A J , where A is given by: 



A = rijA 1 2w n . 



(32) 



Note, that one finds A oc N/tu n , if one measures the 
Zeeman energy of the effective field in units of eni- Ex- 
perimentally, this parameter is related to the ratio of the 
nuclear and the external magnetic fields A and 65, re- 
spectively, giving rise to a small A< 1 for large external 
fields. 

In this parameter regime, all orders higher than second 
order are strongly suppressed and, thus, can be neglected. 
As we show below, neglecting terms of order 0(A) limits 
the range of validity of our analysis to times of order 
t <C A^V///. In order to extent to longer times, it would 
be necessary to include^ higher orders of the self-energy, 
which is however beyond the scope of this article. 

The self-energy exhibits to all orders Cy a 4 x 4 struc- 
ture 



(s 

E?>( S 

eL 2 1(, 



1 

"4 

i 

4 

-4 2) (.s) Ef t ) ( S ) = -4 2) ( S ). 



Tr/ (g^T^g+p! + g-g+^pi) 
Tr/ {g+F^g-pi + g + g-F^pi) 



Tr/ (g+G\g-Pi + g-g+Gif>i) , 
Tr/ (g+Gig+pi + g+g+GtPi) > 
Tr/ {g-G^g-pi + g-g-GiPi) , 
Tr/ (g-Gig+pi + g+g-GtPi) ■ 



(34) 

(35) 
(36) 
(37) 

(38) 

(39) 

(40) 



The super-operators of the electron spin dynamics are 
given by 

F t , l (L + ) = (s-ia u QL+/2)- 1 (41) 



and 



G u (£_) = (s + , (42) 

where 0^4. = ±1. Using these equations in combination 
with Eqs. (22) and (31), we can compute all parts of the 
self-energy for arbitrary anisotropy, which is presented in 
detail in the Appendix. The second order self-energy is 
given by: 



4 2) ( S ) 



= -iN 



-W 



A 2 
A 2 

(43) 



Eg? 00 



A 2 
A 2 

(44) 

S ) = -Sg ) (a), = (45) 



/ E tt (s) E n (s) \ 

=<■> - s f 5) s f ] E+ : w s+ ° (s) ■ (^) =&w - w + w] . 

V E_+(s) E__(s)/ (46) 

sL 2 l(.s) = - W A? -^ X ± [c_7_^ (,s) + c+ / + ^ («)] , 
which shows a block-diagonal form indicating, that the A z 
longitudinal and transverse subspaces are decoupled. In 

(47) 



9 



Z ( +l(a) = iN 



C-I_- x ,(s)+c + I + - X:i (a)] = 



A_A, 
~2Af 



= iJV' 



+ A 
,A_A + 
2AF 
2A_A 



+ 



+ v(2) 



A 2 



(48) 



(49) 



where the coefficients c± arise from the expectation val- 
ues of the operators h± with respect to nuclear mag- 
netic field eigenstates^l. The parameter A z = sign(A 2 ) 
takes into account the effect of the sign of the longitudi- 
nal anisotropy coefficient A z . In the following all calcu- 
lations are executed using A z = 1, since the results for 
X z = — 1 can be easily derived from these computations 
as we show at the end of the following section. 

For a uniformly polarized spin 1/2 system like 
graphene, the coefficients c± take the simple form c± = 
[1 =F p]/2 depending on the polarization p. Moreover, 
the computation of the self energy parts creates super- 
operator matrix elements [£7j] g n with respect to nuclear 
eigenstates which in turn give rise to the functions 



I± (s) = s [log (s=Fi)- log (a)] ± i . 



(50) 



These functions are calculated applying the same contin- 
uum limit used above to obtain Eq. (26|. 

Analyzing the above Eqs. (43) to (49) in more de- 



tail, one finds that all longitudinal self-energy parts have 
always finite values, whereas the off-diagonal transverse 
parts S±ip(s) vanish for A^ = A y . In contrast to this, 
all transverse parts contribute for X x ^ A y , where the 
anisotropy of the HI does not only change the prefac- 
tors, but also gives rise to additional contributions to the 
transverse spin component (S+}(t). 



B. Inverse Laplace transform 

With the specific form of the self-energy in second or- 
der, we are able to calculate the time evolution of the 



electron spin density matrix by solving Eq. ( 28 1 . The ex- 



pectation values for the transverse and longitudinal spin 
components in Laplace space are obtained by 



(51) 



(Sp)(s) = Tr s [Spps(s) 



where /3 = z, ±. The transverse components generally 
consist of two parts 

/- x ( S ±)n &+L (*») (SV) n 

(s± )(* ) = jrf\+ n T ( \„ T\° , (52) 

D±(s) D±(s)D T (s) 
where the denominator functions are given by 



D ± (a) = sTiuJn + isL'i (s) 



(53) 



Symbol 



_L 



(A X , \y, \ Z ) 



(5 + , 2 ) |(5_) ^|(5-) 

(s. tl ) (s z ) (s z ) (s z ) 

(5.,a> 1 I 



TABLE III: Prefactors of the different models for graphene 
in perpendicular and parallel magnetic field and for a general 
anisotropy {A;}. Note that there is no special choice of the 
initial values (S/3 > \) (j of the electron spin. 



With the introduction of the effective nuclear number 
N = N(\1 + Xy)/2Xl, these functions are formally equal 
for arbitrary anisotropy. Hence, the mathematics are 
greatly simplified by the use of these general denomi- 
nators given in Eq. (61) below. Using Eqs. (48) and 



( 49 1 , one can rewrite the transverse spin component in 



terms of two generalized functions: 

(S ± )(s) = (S ±tl )(s) + (S ±>2 )(s) 

_ (S±,l) (£±,2)0 ' ( S±iuJ n) 



D± (a) 



D± (s)D T (a) 



(54) 



where the prefactors for different cases are listed in Tab. 

eh 

In the same manner, we will reformulate the longitu- 
dinal spin component 



<&>(*) = 



D z (a) D z (s) 
where the denominator function reads 



Dz [a) 



.(2) 



» 



E<J (a) 



and the numerator of the second part is given by 



N z (a) = - 



2s 



.(2) 



(a) 



.(2) 



(a) 



(55) 



(56) 



(57) 



Again, we define two more functions 

(S z ) (a) = (S zA ) (a) + (S., 2 ) (a) 



D z (a) 



>D Z (a) 



(58) 



By extending the functional dep ende ncies of all general- 
ized functions given in Eqs. (54) and (58) by uj n and A z , 



one can easily show, that the following relation holds: 

(Sp,i)(s,uj n ,-X z ) = (Sp t i)(s,-u) n ,\ z ) , (59) 
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where ft = ±, z. As a consequence of this, a negative 
coupling A 2 < in the longitudinal direction can be effec- 
tively regarded as an inversion of the external magnetic 
field. 

Moreover, further simplifications are possible by relat- 
ing the longitudinal denominator D z (s) to the denomi- 
nator D + (s) of the transverse spin parts through a shift 
of its variable from s to s — s + iuj n : 



D z (s 



D+ (S) 



N 



s — i2u>n 



0(d) 



(60) 



As we will show below, the relevant part of the complex 
plane is given by |Im(s)| < uj n . Thus, the second contri- 



bution on the right-hand side is nearly a constant of order 
0(A) for all relevant parts of the complex plane. From 
Eq. (61 1 below, it is clear, that this constant is small only 
compared to the denominator D + (s) for |Im(s)| < uj n . 
Thus, this correction to D + (s) will lead to only a slight 
and, hence, irrelevant shift in the following calculations, 
which we will neglect for simplicity. However, reaching 
s « iw n , this approximation, becomes worse, and, hence, 
we will take further efforts to investigate this regime. In 
order to be consistent, we have to shift the numerator 
N z (s) in the same way as the denominator. The shifted 
numerator and the transverse denominator are explicitly 
given by 



D± (s) = s+iuj n ± iN (e_ — c+) 



Ns 



: jln(s-i) - ln(s) J + c±{ In (s + i) - In (s) J 



(61) 



N z (s - kj n ) 



N 



2 s 



i (c_ + c + ) + S 



. { In (s - i) - In (s) } - c+{ In (5 + i) - In (s) } 



+ 0(6) . (62) 



The corrections to N z (s — icj„) are of order O(S) for |s| < 
to n and consequently negligible due to their smallness. In 
the vicinity of s ~ iw n , though, we face similar difficulties 
as for the denominator demanding a careful treatment. 

These peculiarities are especially important for the in- 
verse Laplace transformation, which is applied to recover 
the time-dependence of the electron spin. This transfor- 
mation is obtained by evaluating integrals of the form 



7+100 



<3s,a) (t) 



2?ri 



(S^ x )( S -in p )ds , (63) 



7—100 



along the Bromwich contour by means of complex analy- 
sis, where the indices are given by /? = z, ± and A = 1, 2. 
The frequency shift 



for 13 = ± 
u n for j3 = z 



(64) 



is relevant only for the z component of the electron spin 
and is used to simplify the mathematics as explained 
above. The constant 7 £ M. is chosen, such that all sin- 
gularities have a real part smaller than 7. These singu- 
larities are generated by the denominators D± (s) each 
possessing three zeros at Sj t ± (j — 1,2,3), which come 



in complex-conjugated pairs Sj !+ = 



Furthermore, 



there are three branch cuts from s = 0, ±i to —00, whose 
position in the complex plane is illustrated in Fig. [2] 

In order to evaluate the integral in Eq. (63), one can 
close the contour as depicted in Fig. [2j where the inte- 
gral over the great circle vanishes according to Jordan's 



lemma. Note that the poles with finite real part S2,± are 
outside of this contour and, hence, do not contribute to 



the integral. Therefore, the solution of Eq. (63) gener 



ally consists of residues arising from the remaining poles, 
si,± and S3,±, the integrals Bp x (t) along the upper and 
lower branch-cut as well as an integral Ip \ (t) along the 
imaginary axis: 



[e st -(S^)(s) 



(s P ,x) (t) = £ n< » 

s={si.±,s 3>± } 

1 ■■'"'*{$>!,*(*) + W)}, 



27ri 



where the integrals are explicitly given by 



(65) 



Bl x (t) = lim e iat j dx e xt (Sp,x) {x + ia(l + rf) - iO^) 



(66) 



with a = ±1 and 



Ip* (t) = ij dy ■ (§ PiX ) (iy - iftp) . (67) 
-1 

In order to simplify the notation of the results and to 
make their interpretation easier, we use the short hands 
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Im (.5) 




Re (s) 



FIG. 2: (Color online) Illustration of the contour integral in 
the complex plane with the analytic features of the denom- 
inator functions D+(s) and D-(s) consisting of branch-cuts 
(red dashed lines) and poles (blue and green circles). The 
integral in Eq. ( 63 1 is completed to a closed contour by an 



exponentially vanishing integral over the great circle, integrals 
along the upper and lower branch-cuts, Bp X (t), and an in- 
tegral along the imaginary axis {t). Within this closed 
contour lie the poles Sj,±, j = 1,3, whereas the poles S2,± 
are not encircled and, thus, do not contribute to the integral. 
Note that the poles Sj _ are relevant only for the calculation 
of (&,+>(*). 



N, A, A = A-{Xl + X 2 y )/2X 2 z7 S, and 5 = 5 ■ {\ 2 X + \ 2 y ) / 2 
listed in Table [TTJ as well as the relations c_ + c + = 1 
and c_ — c + = p, which arc fulfilled for nuclei with spin 
I = 1/2 as considered here. 



V. RESULTS 

Inverse Laplace transformation of the first 
transverse spin part 



in 

v 



HI 




FIG. 3: (Color online) Real part of P+.i(t) (solid, black) as 
a function of time obtained by numerical integration. For 
times t > tbi, it is asymptotically described by an oscil- 
lating function (dashed, red), whose amplitude decays with 
~ t^ 1 (dotted, blue), where the amplitude is proportional 
to the polarization p. The oscillations show a frequency of 
/ = tZ\ oc |A z |j4 determined by the HI. 



in 

v 
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ll \ 


t 

sin(Q 
t 














1 ' .** 
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Numerics 
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t [7T-T H |] 



10 



We begin the inverse transformation into the time do- 
main with the first part of the transverse spin component 
(S+.i) (s), which is exemplary for the calculation of all 
other expectation values (Sp t \>{t). 

First, we analyze the residues arising from the poles 
Sj,-|-, which are summarized in Tab. 



IV The pole s 3 . 



located on the imaginary axis gives rise to a purely oscil- 
lating term, where the frequency is given by the effective 
magnetic field u> n . This oscillating part corresponds to a 
simple precession of the electron spin around this mag- 
netic field. 

The calculation of this residue also nicely illustrates in 
which sense the disregard of the fourth order contribution 
£( 4 )(s) oc A of the self-energy sets an upper time limit. 
For simplicity, we assume that this contribution can be 
described by a complex valued constant £( 4 )(s) ps £A 
of order 0(A), where we neglect any dependence on s. 
According to Eq. ( 53 ) , this constant shift can be formally 



treated as a modification of the effective magnetic field 
iu„ — > iuj n — £A giving rise to an additional exponential 
factor exp(— £A-i). Hence, our predictions are valid only 



FIG. 4: (Color online) Imaginary part of P+ t i(t) (solid, black) 
showing a similar behavior as the real part depicted in Fig. [3] 
In contrast to the real part, the amplitude of the imaginary 
part does not depend on the polarization. 



for times t <C /St 1 thi, for which this factor is irrelevant. 
In order to extend this limit, one would have to take the 
fourth order contribution with its full s-dependence into 
account. 

This line of arguing is, however, not directly applica- 
ble for the oscillating portion originating from the pole 
si ; +, since this pole has a more complicated structure. 
It is located near the lower branching point at — i, as is 
illustrated in Fig. [2j This residue generates an ampli- 
tude, which is exponentially small exp[— (c + A) -1 ] <C 1 
in a large magnetic field. As a consequence of its small- 
ness, we will neglect this contribution in the following 
for simplicity. Finally, the pole S2.+ does not contribute, 
because it is outside of the contour. 

Since the integral over the great circle vanishes due to 
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D+ oo = 



Res [e st -(S + ,i)(s)] i 



Res [e st -(S z , 2 )(§-iu J „)] i 



Si,+ : — i 



s 3j+ : ioj„ 



1 + 2 c + e Ac + 



Nc+2 a + 



Mo 



Nc+2 c + 



- e _^_ i( „„ +1)t 



TABLE IV: Zeros of the denominator D + (s) and the corresponding residues for the transverse and longitudinal electron spin 
components (S+,i\(t) and (5*2,2) (i), respectively. The purely imaginary pole 53,+ gives rise to an undamped oscillation around 
the effective magnetic-field uj n . The special case of s = 0, which corresponds to s = S3.+ = iw n , gives rise to a longtime average 
(•Sz)^, which we discuss in more detail in the text. 



Jordan's lemma, the only remaining, unknown expres- 
sions arise from the integrals along the branch-cuts and 
along the imaginary axis. The calculation of these in- 
tegrals is, however, mathematically very challenging due 
to the fact that three different scaling behaviors are in- 
volved. The inverse Laplace transformation itself gives 
rise to an exponential factor exp(si), while the denom- 
inators contain both logarithmic and power law terms 
hampering analytical solutions to these integrals. Nev- 
ertheless, analytical considerations give valuable insights 
to the structure of the results. First, we checked, that 
there are no contributions present in P/3 t \(t), which di- 
verge for longer times and, thus, would lead to unphysical 
results. Moreover, one can testify, that the leading or- 
der contributions arising from the branch cut integrals 
(t) and from the imaginary integrals (t) cancel 
each other leaving terms of order <5 = A 2 /N. While it is 
easy to show that the contributions of this order stem- 
ming from the branch cut integrals are oscillating with 
a frequency r^j determined by the HI, this is not evi- 
dent for the imaginary integral. Furthermore, there is no 
obvious way to analytically extract more information on 
the time dependence of the amplitudes such as the form 
of a possible decay. Hence, we use numerical methods to 
find the results for all integrals, which are subsequently 
summed up in order to give the functions Pp t \(t) defined 
in Eq. ( 65 ) above. 



Eq. (65): 



The function being relevant for (S'+ ! i}(t) is given by 
P+ i(t), whose real and imaginary part is plotted in Figs. 
[3] and [4j respectively. For times t > thi , we find that 
the sum of the branch-cut contributions is asymptoti- 
cally described by an oscillating term, whose amplitude 
is decaying with a power law: 



(S+,i>o 



= -27Ti(5 



sin(i) . cos(i) 



vp- 



(68) 



The final result for (S'+ ! i)(t) is obtained by summing 
up the residues and power law contributions according to 



(S+,i) (*) 



1 + (A 2 + A 2 ) 



S 



2 



sin(t) 
t 



lp- 



cos(t) 



(69) 



where we reintroduced the explicit dependence on the 
anisotropy of the HI using the relations summarized in 
Tab [TTJ Note, that a negative coupling X z < in z- 
direction can be easily handled using Eq. ( 59 ) . 



Fo r th e special case of isotropy within the x-y plane, 
Eq. (54) readily gives (S+, 2) = 0. Thus, setting X x = X y 
in the above equation already allows us to interpret the 
dynamics of this type of systems, which are physically re- 
alized in GaAs or graphene subjected to a perpendicular 
magnetic field. According to Eq. ( 69 ) only a small frac- 
tion of order 6 = A 2 /N of the initial transverse spin de- 
cays and does so in a power law, while most of the trans- 
verse spin is preserved and oscillates with a frequency 
determined by the effective magnetic field uj n . Physi- 
cally this means, that the transverse spin is only little 
affected by the virtual spin flip-flops, which are strongly 
suppressed because of the enormous energy difference be- 
tween the Zeeman splitting and the HI energy. 



B. Inverse Laplace transformation of the second 
transverse spin part 

So far, we have calculated the time-dependence of the 
first transverse spin part (S'+.i)(i), which fully describes 
the behavior of a system with rotational symmetry in the 
transverse x-y-plane. For a system having a broken ro- 
tational symmetry with A^ =/= X y , one has additionally to 
calculate the second transverse spin function (S , +j2 )(t), 
which contains two denominator functions D+(s) and 
D_(s) instead of one. 

Turning first to the residues, we can take advantage 
of the fact, that the zeros of these denominators come 
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FIG. 5: (Color online) Real part of P+^ii) (solid, black) as 
a function of time obtained by numerical integration. For 
times t > thi, it is asymptotically described by an oscillating 
power law decay (dashed, red; dotted, blue) ~ t~ , where 
the amplitude is proportional to the polarization p and the 
frequency is given by / = t^) ■ 



in complex conjugated pairs = s* _ because of the 
relation D*As) — D-(s*). Therefore, the calculation of 
the residues is straightforward, resulting in the contri- 
butions listed in Tab. |Vj The poles Si t ± generate expo- 
nentially suppressed terms, which are equivalent to the 
si,+-residue of (5 +i i)(t). This means in particular, that 
these residues are also negligible for large magnetic fields. 
In contrast to these paired terms, the pole S3 t + creates 
a purely oscillating component, while the residue of pole 
vanishes identically. 
The calculation of the sum of the branch-cut integrals, 
P+^i), i s again obtained by numerical integration ac- 
companied by analytical considerations. We find, that 
only the real part of P +i 2(i) has a relevant contribution 
of order 0(6), while the imaginary part is much smaller 
of order 0(S/ui n ) and, thus, neglected. The time evolu- 
tion of this real part is, up to a factor of two, analogous 
to the real part of P + i(t), which can be read off from 
Fig. [5} 

We find, that the branch-cut contribution can be 
asymptotically described by 



P+,2 (*) 



= -2iri5 



- \2p 



cos(i) 



(70) 



for times t > thi- In order to obtain the full (S + ^)(t) 
term, we sum up the residues and power law contribu- 
tions yielding: 



(S + ,2)(t) 



Auj n t 



\2p 



cos(i) 



t 



(71) 



In combination with the result for (S+.i)(i) given in 
Eq. (69), we are now able to formulate the time depen- 



dence of the transverse electron spin component (S+)(t) 



for arbitrary anisotropy and a general initial condition 

(S±) = (S x ) ±i(Sy) Q : 

(S+)(t) = (S + , 1 )(t) + (S + , 2 )(t) 

= ((s x ) + i(s y ) ) [l + (\l + Xl) S - 



+ (Xl(S x ) + iXl(S y ) Q )5 i 

-^l(s x ) Q + ixl(s y ) )s P c -^ 



(72) 



In the limit X y — X x , this result reproduces the previ- 



ous result given in Eq. ( 69 ) . As in the isotropic limit 



discussed above, we find that due to the large Zeeman 
splitting most of the transverse electron spin is preserved 
and precesses around the effective magnetic field uj n . The 
decaying part, however, differs from the isotropic case 
and allows us to analyze how the effect of the HI changes 
as X x =^ X y . As exp ecte d for the transverse expectation 
value (S* + )(t), Eq. (I72J) is formally invariant under the 
exchange x •<-> y, which allows us to discuss the result for 
a specific choice of > X y without loss of generality. 
Modifying the coupling until X y /X x = is reached, the 
amplitude of the sinusoidal part is more and more domi- 
nated by the initial electron spin in y direction, while the 
amplitude of the cosine term is increasingly governed by 
the initial component (S X } Q . The reason for this is the 
broken rotational symmetry in the x-y-plane due to the 
anisotropy of the HI. This in turn makes the time evolu- 
tion of the electron spin dependent on its initial prepa- 
ration. 

Furthermore, the amplitude of the cosine term also de- 
pends on the polarization p of the nuclear bath, which 
quantifies the excess of one nuclear spin orientation 
over the other. For increasing p one type of the HI- 
induced scattering processes, for instance h + , becomes 
more likely, while the other one is suppressed, since its 
phase space is more and more limited. Presumably, this 
can explain the polarization dependence of our result. 

However, within the Nakajima-Zwanzig formalism, it is 
not possible to single out microscopic processes explain- 
ing the specific form of the dependence on the anisotropy 
and the polarization demanding for successive studies us- 
ing different techniques. 



C. Longitudinal electron spin for arbitrary HI 

Next, we turn to the inverse Laplace transformation 
of the longitudinal spin components, where we use the 
shifted denominator and numerator functions given in 
Eqs. (60) and (62), respectively, to re-express Eq. (63) 



in terms of the new coordinate s: 



(S,,l/2) (t) 



-\UJ n t 



7+ioo 



27ri 



{S x , 1/2) (s-iw„)ds 



(73) 
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D± (s) = 



Res [c 3t (S + , 2 )(s)] s=Sjt+ 



Res [e st <5- + ,2>( S )]^ 



S — S t _ — S ■ 



Si, 



S3, 



<s + , 2 ) [1 + 



"( 5 '2,+ ) 



ATc+ 2 c + 



TABLE V: Residues of the second transverse spin part (S2, +)(£)• The poles Sj,+ are listed in Tab. IV Here the ss^-poles 
produce inequivalent residues due to the rewritten form of (S , 2,+) (s) in Eq. (541. 
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FIG. 6: (Color online) Real part of P Z} i{t) (solid, black) as 
a function of time obtained by numerical integration showing 
a similar behavior as Re[P+,i(t)] depicted in Fig. [3] In con- 
trast to this, however, its amplitude does not depend on the 
polarization. 



FIG. 7: (Color online) Imaginary part of P z ,2(t) (solid, black) 
as a function of time obtained by numerical integration show- 
ing a similar behavior as Im[P+ I i(t)] depicted in Fig. [4] De- 
viating from this, the amplitude depends on the polarization 
of the nuclear bath. 



With this reformulation, the integrand (S Zt -u(s — iu> n ) 
has the same analytic structure like (S+,i)(s) for 
Im(s) <C u) n , and, hence, we can take advantage of pre- 
vious results. Therefore, the contributions of the poles 
si j+ as well as the branch-cut and imaginary integrals 
are readily obtained by replacing the prefactor (iS i ;+ ) n 

by (S z ,i) ■ exp[-iw n i] in Tab. [iv|and in Eq. (|6S|) . As 
we discussed extensively in context wi 
Eqs 



can be described by the following function: 



with the derivation of 
(|60|) and (|62|), the calculation of the residue of the 



pole S3.+ = \u> n needs a more sophisticated treatment, 
which will be done in detail below. 

In order to find the full time dependence of the longi- 
tudinal spin component, we also need to compute the 
second contribution \S Z ,2) (i), which follows the same 
procedure as (5 + .i)(t). This calculation is in principle 
easy, but lengthy and is, therefore, not presented in de- 
tail. The residue of the pole si j+ is listed in Tab. IV 
whereas the outcome of the pole at S3 + = iuj n will be 
also presented in more detail below. 

The integrals making up P z ,2(t) are again calculated 
numerically. The time dependence of these integrals is 
presented in Figs. [6] and [7j from which their power law 
decaying behavior becomes obvious. Asymptotically this 



P,,a(t) S 

= Z7T1 — 

2 



sin(t) cos(t) 
p— 1 — : — 



(74) 



Note, that this result is smaller by a factor of two com- 
pared to the outcome of P +: i(t) given in Eq. (68) 



and that the real part depends on the polarization in- 
stead of the imaginary part. At last, we have to evalu- 
ate the residues at the pole £3,+ for both (S Zt i)(t) and 
(S Zt 2)(t). If we shift back to the original coordinate sys- 
tem s — > s = s — VjJ n , this residue can be rewritten as 
a limit for long times according to the properties of the 
Laplace transform: 

Res ^-'^{S^is-iun) 
= Urn s e st (S zA/2 ) (s) = (S z , 1/2 ) . (75) 

This longtime average of the z-spin component \S Z ) — 
(5*2,1) + (S z ,2) is actually calculated by performing a 
series expansion in (s±iw n ) 1 of Eq. (58) and subse- 
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quently taking the limit s 0, which finally leads to: 

(U+if(U+ (5) 



QD with an in-plane magnetic field is exemplary for the 



1 + S + O 

(S z ) + X x X y ^ 
1 + ^S 



N 



(76) 



With this result, all contributions to the longitudinal 
electron spin component are known: 



{S I )(() = (S,)+^e- 1 ™-'- 

>,.)„ -§<*,>.} 

-■{»<4..>. + j<M„}2*i 

\2 I \2 



Air. A 



sin(t) 



A, 2 , J * 



(77) 



which consists of a constant contribution (S z ) and a 
small decaying portion of order 0(8). 

In contrast to the transverse spin part, the power law 
decaying term oscillates with two frequencies uj n , deter- 
mined by the effective magnetic field, and t^\, which is 
set by the strength of the HI. The amplitude of this os- 
cillation is decaying similarly to the corresponding part 
of (S+)(i) with a power law ~ t -1 . As a consequence 
of this, the longitudinal spin component settles to a con- 
stant value \S Z ) for longer times, which is up to small 
corrections of order 6 given by the initial value. 

Note that all results are valid for arbitrary coupling 
constants A^, X y and A z , which particularly implicates, 
that these findings are true for a graphene QD with both 
a perpendicular and a parallel magnetic field. A possible 
negative sign of the longitudinal coupling constant X z < 
0, as it is the case for a parallel field, can be handled by 
changing the sign of the magnetic fields, while using |A Z |, 
which is allowed by symmetries of (S z ,i/2){t) described 
inEq. |59|. 



D. Discussion of the results 

Motivated by the physics of the HI in a graphene QD, 
we have analyzed the effect of an anisotropic HI, where we 
find two distinct classes, which are characterized by A^ = 
A y and A^ ^ X y , respectively. Representative for the for- 
mer is a graphene QD subjected to a perpendicular ori- 
ented magnetic field, where — 2X X = —2X y = X z = 1, but 
also the earlier investigated^ isotropic model for GaAs is 
a member of this class. In contrast to this, a graphene 



-2Ai, 



-2X, 



1. 



second class of HI, where X x 

From Eqs. (72) and (77 1, one sees, that the overall dy- 



namics of the electron spin is the same for both classes of 
anisotropy. The transverse components exhibit a dom- 
inant oscillating contribution, which describes a simple 
precession of the electron spin around the effective mag- 
netic field consisting of the external and nuclear fields. 
Moreover, the transverse component features another 
term arising from the Hl-induced flip-flop processes. Due 
to the large Zeeman splitting of the electron spin states, 
these processes are suppressed leading to a very small 
absolute value of this contribution compared to the pre- 
cession term. Additionally, its amplitude oscillates with 
a much smaller frequency set by the HI and decays as a 
function of time, which is asymptotically well described 
by a power-law for times t > Thi, while for even shorter 
times only numerical results are within reach. The lon- 
gitudinal spin component exhibits a similar power law 
decay, whose amplitude is in addition modulated by an 
oscillation with the frequency corresponding to the effec- 
tive magnetic field. The main part of the longitudinal 
spin is, however, preserved and is given up to small cor- 
rections by the initial value of the electron spin. 

While these general remarks on the electron spin dy- 
namics are true for arbitrary anisotropy, there are, how- 
ever, both qualitative and quantitative ramifications of 
different choices of the couplings Xi , which is reflected in 
the amplitudes of the decaying contributions. According 
to Eq. (771 the longitudinal part shows only a quanti- 



tative dependence, whereas the transverse spin parts are 
affected quantitatively as well as qualitatively by different 
types of anisotropy as can be seen in Eq. (72). For the 



first class of anisotropy characterized by X x = X y , we find 
that all amplitudes are proportional to the initial value 
(S+) = (^x) + 'i(Sy) Q leading to pure quantitative dis- 
tinctions within this class. Since the qualitative concur- 
rence arises from A^ = Aj, independently of the specific 
values of the constants, it seems that rather generally, an 
anisotropy between the transverse and longitudinal spin- 
subspace, does not lead to different physical behavior. In 
contrast to this, the second class with A^ ^ Aj, exhibits 
a X x and Aj, weighted mixing of the initial values \S x / 
and \Sy) . This is caused by the broken rotational sym- 
metry the x-y plane defining a precise reference frame 
according to which \S X J Q and (S y ) can be measured. 

Despite this mixing of amplitudes, one can addition- 
ally expect different timescales of the HI as is discussed 
in Sec. [XTJ. For the parallel and perpendicular case, one 
finds Tjjj/tjjj = 2 due to the different coupling X z in 
the longitudinal direction. Thus, longer coherence times 
should be measurable for an in-plane magnetic field. 



VI. CONCLUSION AND OUTLOOK 

In this article, we have shown, how the analysis of the 
HI by means of a generalized master equation can be ex- 



16 



tended from the isotropic case^to a general HI, which can 
in principle be completely anisotropic (X x 7^ X y 7^ X z ). 
Investigating a graphene QD, with an external magnetic 
field applied either perpendicular or parallel to the car- 
bon plane, we were able to study two specific realiza- 
tions of an anisotropic HI, which are representative for 
two different classes of anisotropy. While the former sys- 
tem describes a situation, where the anisotropy exists 
between the longitudinal and transverse spin subspaces 
(\ x = X y 7^ A z ), the latter exhibits the anisotropy within 
the transverse subspace (X x 7^ X y ), where no other spe- 
cial relation between X X ,X V , and X z is assumed. As we 
showed above, the first class of anisotropy gives rise only 
to quantitative differences leaving the analytical struc- 
ture of the electron spin dynamics unchanged with re- 
spect to the isotropic case. In this model, most of the 
spin amplitude is preserved, while a small portion decays 
with a non-exponential behavior. Presumably, the reason 
for this universality is the fact, that the rotational sym- 
metry within the transverse subspace is not broken for 
X x = X y . This symmetry is, however, not preserved for 
the second class of anisotropy, in which 7^ X y , leading 
to a different amplitude of the power law decaying con- 
tribution. While its order of magnitude is unchanged, we 
find an anisotropy weighted mixing of the initial ampli- 
tude of the electron spin. As a consequence of this, its 
preparation with respect to the precisely defined x and 
y direction matters in contrast to a rotational symmetric 
system. 

All these results were obtained for times t <C A~ 1 thj, 
which are determined by the ratio of the nuclear and 
external magnetic fields, A cx A/hrfsB z , In this range 
of time, the electron nuclear spin system is in a non- 
Mar kovian regime. Considering even longer times t 3> 
A _1 T/f/, which are not captured by our (second order) 
treatment, the system can, however, return to a Marko- 



vian regime again, as was founcti^U for an electron spin 
in a GaAs QD. Thus, it should be interesting in the fu- 
ture to study how the time evolution of a spin system 
with anisotropic HI behaves for longer times. 

Our findings are qualitatively valid for other systems 
fulfilling the requirements of our model, where the most 
important demands are a Gaussian-like envelope func- 
tion, slow dynamics of the nuclear bath and a sufficiently 
large Zeeman-splitting with respect to the HI energy 
scale. Graphene-specific quantitative changes of our re- 
sults arise from the small HI coupling constant A and 
the low natural abundance nj ps 0.01 of 13 C leading 
to a prolonged time scale of the HI and a reduced nu- 
clear magnetic field, which make the design of qubits less 
challenging, since computation-cycles can last for longer 
times and lower external fields are sufficient. These latter 
effects would become even more important in isotopically 
purified samples featuring nj < 0.01. However, if one in- 
creasingly reduces the amount of 13 C in the graphene 
QD, one will reach a regime, where it is not appropri- 
ate anymore to think of a bath of nuclear spins. Unfor- 
tunately, this regime is not accessible by our treatment 
calling for successive studies, which deal with a finite 
number of nuclear spins^H where in particular a more 
detailed study of the influence of the initial conditions 
should be possible. 



Acknowledgement 

We acknowledge support from the Priority Program 
1459 "Graphene" of the DFG and from the Euro- 
GRAPHENE Program of the ESF. Furthermore we 
would like to thank Bill Coish, Patrik Recher, and Stefan 
Walter for valuable discussions. 



1 D. Loss and D. P. DiVincenzo, Phys. Rev. A 57, 120 (1998) 

2 V. Cerletti, W. A. Coish, O. Gywat, and D. Loss, Nan- 
otechnology 16, R27 (2005) 

3 J. Schliemann, A. V. Khaetskii, and D. Loss, J. Phys. Con- 
dens. Matter 15, R1809 (2003) 

4 W. A. Coish and J. Baugh, Phys. Status Solidi (B) 246, 
2203 (2009) 

5 L. Cywinski, Acta Phys. Pol. A 119, 576 (2011) 

6 J. Schliemann, A. V. Khaetskii, and D. Loss, Phys. Rev. 
B 66, 245303 (2002) 

7 A. V. Khaetskii, D. Loss, and L. Glazman, Phys. Rev. Lett. 
88, 186802 (2002) 

8 A. V. Khaetskii, D. Loss, and L. Glazman, Phys. Rev. B 
67, 195329 (2003) 

9 W. A. Coish and D. Loss, Phys. Rev. B 70, 195340 (2004) 

D. Klauser, W. A. Coish, and D. Loss, Phys. Rev. B 73, 
205302 (2006) 

1 D. Stepanenko, G. Burkard, G. Giedke, and A. Imamoglu, 
Phys. Rev. Lett. 96, 136401 (2006) 

2 G. Giedke, J. M. Taylor, D. DAlessandro, M. D. Lukin, 



and A. Imamoglu, Phys. Rev. A 74, 032316 (2006) 

13 W. Yao, R.-B. Liu, and L. J. Sham, Phys. Rev. B 74, 
195301 (2006) 

14 D. Klauser, W. A. Coish, and D. Loss, Phys. Rev. B 78, 
205301 (2008) 

15 W. A. Coish, J. Fischer, and D. Loss, Phys. Rev. B 77, 
125329 (2008) 

16 J. Fischer, W. A. Coish, D. V. Bulaev, and D. Loss, Phys. 
Rev. B 78, 155329 (2008) 

17 C. Deng and X. Hu, Phys. Rev. B 78, 245301 (2008) 

18 L. Cywinski, W. M. Witzel, and S. Das Sarma, Phys. Rev. 
Lett. 102, 057601 (2009) 

19 C. Testelin, F. Bernardot, B. Eble, and M. Chamarro, 
Phys. Rev. B 79, 195440 (2009) 

20 J. Fischer and D. Loss, Phys. Rev. Lett. 105, 266603 (2010) 

21 W. A. Coish, J. Fischer, and D. Loss, Phys. Rev. B 81, 
165315 (2010) 

22 E. Barnes, L. Cywinski, and S. Das Sarma, Phys. Rev. B 
84, 155315 (2011) 

23 R. Hanson, L. P. Kouwenhoven, J. R. Petta, S. Tarucha, 



17 



and L. M. K. Vandersypen, Rev. Mod. Phys. 79, 1217 
(2007) 

24 D. J. Reilly, J. M. Taylor, J. R. Petta, C. M. Marcus, M. P. 
Hanson, and A. C. Gossard, Science 321, 817 (2008) 

25 X. Xu, W. Yao, B. Sun, D. G. Steel, A. S. Bracker, D. Gam- 
mon, and L. J. Sham, Nature 459, 1105 (2009) 

26 C. Latta, A. Hogele, Y. Zhao, A. N. Vamivakas, 
P. Maletinsky, M. Kroner, J. Dreiser, I. Carusotto, 
A. Badolato, D. Schuh, W. Wegscheider, M. Atature, and 

A. Imamoglu, Nat. Phys. 5, 758 (2009) 

27 I. T. Vink, K. C. Nowack, F. H. L. Koppens, J. Danon, 
Y. V. Nazarov, and L. M. K. Vandersypen, Nat. Phys. 5, 
764 (2009) 

28 M. Gullans, J. J. Krich, J. M. Taylor, H. Bluhm, B. I. 
Halperin, C. M. Marcus, M. Stopa, A. Yacoby, and M. D. 
Lukin, Phys. Rev. Lett. 104, 226807 (2010) 

29 M. Issler, E. M. Kessler, G. Giedke, S. Yelin, I. Cirac, 
M. D. Lukin, and A. Imamoglu, Phys. Rev. Lett. 105, 
267202 (2010) 

30 E. Togan, Y. Chu, A. Imamoglu, and M. D. Lukin, Nature 
478, 497 (2011) 

31 M. S. Rudner, L. M. K. Vandersypen, V. Vuletic, and L. S. 
Levitov, Phys. Rev. Lett. 107, 206806 (2011) 

32 S. G. Carter, S. E. Economou, A. Shabaev, and A. S. 
Bracker, Phys. Rev. B 83, 115325 (2011) 

33 A. Greilich, S. E. Economou, S. Spatzek, D. R. Yakovlev, 

D. Reuter, A. D. Wieck, T. L. Reinecke, and M. Bayer, 
Nat. Phys. 5, 262 (2009) 

34 B. Eble, C. Testelin, P. Desfonds, F. Bernardot, A. Baloc- 
chi, T. Amand, A. Miard, A. Lemaitre, X. Marie, and 
M. Chamarro, Phys. Rev. Lett. 102, 146601 (2009) 

35 S. Foletti, H. Bluhm, D. Mahalu, V. Umansky, and A. Ya- 
coby, Nat. Phys. 5, 903 (2009) 

36 T. D. Ladd, D. Press, K. DeGreve, P. L. McMahon, 

B. Friess, C. Schneider, M. Kamp, S. Honing, A. Forchel, 
and Y. Yamamoto, Phys. Rev. Lett. 105, 107401 (2010) 

37 D. Press, K. D. Greve, P. L. McMahon, T. D. Ladd, 

B. Friess, C. Schneider, M. Kamp, S. Honing, A. Forchel, 
and Y. Yamamoto, Nature Photonics 4, 367 (2010) 

38 H. Bluhm, S. Foletti, I. Neder, M. S. Rudner, D. Mahalu, 
V. Umansky, and A. Yacoby, Nat. Phys. 7, 109 (2010) 

39 H. Bluhm, S. Foletti, D. Mahalu, V. Umansky, and A. Ya- 
coby, Phys. Rev. Lett. 105, 216803 (2010) 

40 C. Stampfer, J. Giittinger, F. Molitor, D. Graf, T. Ihn, 
and K. Ensslin, Appl. Phys. Lett. 92, 012102 (2008) 

41 L. A. Ponomarenko, F. Schedin, M. I. Katsnelson, R. Yang, 

E. W. Hill, K. S. Novoselov, and A. K. Geim, Science 320, 
356 (2008) 

42 S. Schnez, F. Molitor, C. Stampfer, J. Giittinger, I. Sho- 
rubalko, T. Ihn, and K. Ensslin, Appl. Phys. Lett. 94, 
012107 (2009) 

43 F. Molitor, S. Droscher, J. Giittinger, A. Jacobsen, 

C. Stampfer, T. Ihn, and K. Ensslin, Appl. Phys. Lett. 
94, 222107 (2009) 

44 J. Giittinger, C. Stampfer, T. Frey, T. Ihn, and K. Ensslin, 
Phys. Status Solidi (B) 246, 2553 (2009) 

45 X. Liu, J. B. Oostinga, A. F. Morpurgo, and L. M. K. 
Vandersypen, Phys. Rev. B 80, 121407 (2009) 

46 F. Molitor, H. Knowles, S. Droscher, U. Gasser, T. Choi, 
P. Roulleau, J. Giittinger, A. Jacobsen, C. Stampfer, 
K. Ensslin, and T. Ihn, Europhys. Lett. 89, 67005 (2010) 

47 L.-J. Wang, G. Cao, T. Tu, H.-O. Li, C. Zhou, X.-J. Hao, 
Z. Su, G.-C. Guo, H.-W. Jiang, and G.-P. Guo, Appl. Phys. 
Lett. 97, 262113 (2010) 



48 X. L. Liu, D. Hug, and L. M. K. Vandersypen, Nano Lett. 
10, 1623 (2010) 

49 J. Giittinger, J. Seif, C. Stampfer, A. Capelli, K. Ensslin, 
and T. Ihn, Phys. Rev. B 83, 165445 (2011) 

50 S. Fringes, C. Volk, C. Norda, B. Terres, J. Dauber, S. En- 
gels, S. Trellenkamp, and C. Stampfer, Phys. Status Solidi 
(B) 248, 2684 (2011) 

51 P. Recher and B. Trauzettel, Nanotechnology 21, 302001 
(2010) 

52 H. O. H. Churchill, A. J. Bestwick, J. W. Harlow, F. Kuem- 
meth, D. Marcos, C. H. Stwertka, S. K. Watson, and C. M. 
Marcus, Nat. Phys. 5, 321 (2009) 

53 S. J. Chorley, G. Giavaras, J. Wabnig, G. A. C. Jones, 
C. G. Smith, G. A. D. Briggs, and M. R. Buitelaar, Phys. 
Rev. Lett. 106, 206801 (2011) 

54 T.-F. Fang, Q.-f. Sun, and H.-G. Luo, Phys. Rev. B 84, 
155417 (2011) 

55 A. Eichler, M. Weiss, and C. Schonenberger, Nanotechnol- 
ogy 22, 265204 (2011) 

56 A. Morello, J. J. Pla, F. A. Zwanenburg, K. W. Chan, K. Y. 
Tan, H. Huebl, M. Mottonen, C. D. Nugroho, C. Yang, 
J. A. van Donkelaar, A. D. C. Alves, D. N. Jamieson, C. C. 
Escott, L. C. L. Hollenberg, R. G. Clark, and A. S. Dzurak, 
Nature 467, 687 (2010) 

57 J. J. L. Morton, D. R. McCamey, M. A. Eriksson, and S. A. 
Lyon, Nature 479, 345 (2011) 

58 O. V. Yazyev, Nano Lett. 8, 1011 (2008) 

59 J. Fischer, B. Trauzettel, and D. Loss, Phys. Rev. B 80, 
155401 (2009) 

60 The order of magnitude of the HI is also in agreement 
with the results obtained by for graphene flakes 
and an estimation^ based on measurements of the chem- 
ical shift. There are, however, measurements on carbon 
nanotubes, which claim to see a much larger constant of 
A « 100 ne\^2l. 

61 B. Dora and F. Simon, Phys. Status Solidi (B) 247, 2935 
(2010) 

62 J. L. Ivey and R. L. Mieher, Phys. Rev. B 11, 849 (1975) 

63 L. V. C. Assali, H. M. Petrilli, R. B. Capaz, B. Koiller, 
X. Hu, and S. Das Sarma, Phys. Rev. B 83, 165301 (2011) 

64 D. Brunner, B. D. Gerardot, P. A. Dalgarno, G. Wiist, 
K. Karrai, N. G. Stoltz, P. M. Petroff, and R. J. Warbur- 
ton, Science 325, 70 (2009) 

65 S. Saikin and L. Fedichkin, Phys. Rev. B 67, 161302 (2003) 

66 W. M. Witzel, X. Hu, and S. Das Sarma, Phys. Rev. B 76, 
035212 (2007) 

67 T. Ohta, A. Bostwick, T. Seyller, K. Horn, and E. Roten- 
berg, Science 313, 951 (2006) 

68 E. V. Castro, K. S. Novoselov, S. V. Morozov, N. M. R. 
Peres, J. M. B. Lopes dos Santos, J. Nilsson, F. Guinea, 
A. K. Geim, and A. H. Castro Neto, Phys. Rev. Lett. 99, 
216802 (2007) 

69 Y. Zhang, T.-T. Tang, C. Girit, Z. Hao, M. C. Martin, 
A. Zettl, M. F. Crommie, Y. R. Shen, and F. Wang, Nature 
459, 820 (2009) 

70 P. Recher, J. Nilsson, G. Burkard, and B. Trauzettel, Phys. 
Rev. B 79, 085407 (2009) 

71 B. N. Szafranek, D. Schall, M. Otto, D. Neumaier, and 
H. Kurz, Appl. Phys. Lett. 96, 112103 (2010) 

72 E. Fick and Giinter Sauermann, The Quantum Statis- 
tics of Dynamic Processes, 1st ed. ( Springer- Verlag Berlin- 
Heidelberg, 1990) 

73 H.-P. Breuer and F. Petruccione, The Theory of Open 
Quantum Systems (Oxford University Press, Oxford, 2002) 



18 



P. G. Silvestrov and K. B. Efetov, Phys. Rev. Lett. 98, 
016802 (2007) 

A. De Martino, L. DelPAnna, and R. Egger, Phys. Rev. 
Lett. 98, 066802 (2007) 

B. Trauzettel, D. V. Bulaev, D. Loss, and G. Burkard, Nat. 
Phys. 3, 192 (2007) 

A. Matulis and F. M. Peeters, Phys. Rev. B 77, 115423 
(2008) 

M. Titov, P. M. Ostrovsky, I. V. Gornyi, A. Schuessler, 
and A. D. Mirlin, Phys. Rev. Lett. 104, 076802 (2010) 
W. F. Banholzer and T. R. Anthony, Thin Solid Films 212, 
1 (1992) 

F. Simon, C. Kramberger, R. Pfeiffer, H. Kuzmany, 
V. Zolyomi, J. Kiirti, P. M. Singer, and H. Alloul, Phys. 
Rev. Lett. 95, 017401 (2005) 

A. M. Stoneham, Theory of Defects in Solids (Oxford Uni- 
versity Press, New York, 1985) 

A. Abragam, The Principles of Nuclear Magnetism (Ox- 
ford University Press, New York, 1989) 

C. Goze-Bac, S. Latil, P. Lauginie, V. Jourdain, J. Conard, 
L. Duclaux, A. Rubio, and P. Bernier, Carbon 40, 1825 
(2002) 

A. Palyi and G. Burkard, Phys. Rev. B 80, 201404 (2009) 
C. Ramanathan, Appl. Magn. Reson. 34, 409 (2008) 



for arbitrary operators 0%^ and a function f[C±], which 
can be expanded in powers of (anti-)commutators £_ 
(£+)■ Therefore, the result for (£+) is independent 
from the position, at which these super-operators are 
placed, while moving G-\,\. (•£_) to a neighboring position 
involves a change of sign of its argument £_ — > — C-. 
As a consequence of this, it is sufficient to calculate the 
following traces ^i(J~\,]h±h^px) and r £ri{Q^^h±h^pi), 
which give rise to sums over expectation values with re- 
spect to nuclear eigenstates \q) = fe |m£): 

Trj U 7 i h±h T p I S J 
= £<Pl^{l«) <9l^l r > (r\}K\n) (n\p) 

= X>l*i{lff) (r\}(q\h ± \r) h T \n) 



q,r 



q,r 



(q\ h ± \r) (n\q) (r\ h T \n) 



(A3) 



Appendix A: Calculation of the self-energy in 
second order 



In this appendix we discuss the computat ion of t he s elf- 
energy matrix-elements given in Eqs. (34) to (40) in 
more detail. We assume, that the nuclear magnetic field 
operators g± for a specific form of the anisotropy are 
already inserted. Thus, all parts of the self-energy are 
linear combinations with A^-dependent prefactors, where 
the summands contain one super-operator J-^,± or G-f,],, 
two bare nuclear magnetic field operators h± and the 
nuclear initial state pj = \n) (n\. Calculating expectation 
values with respect to the nuclear state \n), all summands 
featuring squared operators h± vanish identically, which 



reduces the number of contributions. In the following, 
we will present, how the remaining expressions can be 
evaluated in a general approach. For simplicity, we will 
neglect the prefactors in the following presentation. 

Due to the linearity and the cyclicity of the trace, all 
self-energy parts in second order can be written as a lin- 
ear combination of terms of the form 



and analogously for the super-operators Gi- Using Eqs. 
( 41 ) and ( 42 ) , one can calculate the expectation values 



of these super-operators: 
{ s - iat,4.^ 



Y 4^ 



A k [ml + ml) ] } 



{ s + ia u /3^£A fc (m£-mD} \ 



(A4) 



(A5) 



where ctf t ± = ±1 and (3 = ±1. The factor A z is the 
anisotropy coefficient in z-direction. Next, we calculate 
the magnetic field expectation values, where we begin 
with the action of a local operator 1^ ± at an arbitrary 
site k in order to simplify later steps: 



fe,± ( 



K) - J 1(1 + 1) - ml{ml ± 1) \m\ ± 1) 



Tr/ (---pi) 



(Al) 



where the free spots can be filled with one raising opera- 
tor h.\-, one lowering operator /i_ and one of the four 
super-operators J^^(£ + ) and G^,i(£-) given in Eqs. 
(41) and (42), respectively. Altogether, this leads to 
3 • 2 • 4 = 24 possible combinations. This number can 
be reduced by using the relational of Liouvillian-like op- 
erators: 



Tr 



(/[A. 



±]Oi0 2 ) =Tc x (0 1 f[±C±]0 3 



(A2) 



M ± K) 



K±l)(g)K), (A6) 

l^k 



We introduced the shorthand notation M±(m1), which 
obeys the relation M±(m q k =Fl) = M^(m q k ). Note, that by 
the action of Ik,± only one single local state was changed 
while all other local states remain unchanged. Using the 
above equation, one can calculate the expectation value 
of the nuclear magnetic field operators h± with respect 
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to eigenstates \p), \q) of the /i z -component: 



(P\ h± \q) -- 

= (p\^2A k I k , ± <g>\ml) 
k p 

= $>* M ±K) (pi { K ± i) ® K) } 

= E ^ M ± K) { (8) K I } { K ± 1) ® (g) |mf) } 

k V tjtk 

= E A k M±(m 9 k ) 5 m », m »±i 5mf„mf< ( A7 ) 

He Z,i'#fc 

This equation totally sets the relation between the two 

r / l N 

sets of product states < 9 > . Inserting this result in 
Eq. dAl, we find: 



r 

= {s + ia n /3y E ( TO fei - m fei) } 
fci 

&2 l.l'^ki 



y AlM*(ml) 
s±ia n /3^A k2 



(A8) 



The remaining two equations are obtained in the same 
manner leading to: 



E^ w 

r 

E 



hi 



hi 



s - iat,4./3 (w„ ± 4f A fc2 ) ' 



(A9) 



where we used, that 
A 



j E ( m k' + m k>) n ^r.K^K.™^ 1 



l^k 



\ z (h z ) nT ^A k 



(A10) 



The functions M\(m^) over nuclear eigenvalues can be 
replaced by their average assuming a nuclear state which 
is highly degenerat e! 9 * 22 ^: 



«M 2 (m£)))= c± . 



(All) 



Finally, for a large nuclear spin system with N 3> 1, 
the remaining sums in Eqs. (A8) and (A9) can be re- 



placed by integrals in the continuum limit. Changing 
to dimensionless units by measuring energies in units of 
(hi = \\ z \niA/2N , one finds up to small corrections^. 



Trj 
4N 



(jf4 (/?£+) h±h T pi^ 

,2 C T / T {a T ,/3A.}( S " i«t4^n) (A12) 



(A13) 



4N 



where A z = A Z /|A Z | and 

i± (s) = s [log (s T i) - log 0)] ± i . 



(A14) 



Applying this continuum limit, however, sets an upper 
bound t -C <J N/2thi as discussed in the main text. 
Knowing the four basic terms given in Eqs. (A12| and 



(A13), respectively, all other remaining possible sum- 



mands to the self-energy are readily obtained by using 
Eq. |A2|. 



